Even-Odd Correlation Functions on an Optical Lattice 
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We study how different many body states appear in a quantum gas microscope, such as the one 
developed at Harvard [Bakr et al. Nature 462, 74 (2009)], where the site-resolved parity of the atom 
number is imaged. We calculate the spatial correlations of the microscope images, corresponding 
to the correlation function of the parity of the number of atoms at each site. We produce analytic 
results for a number of well-known models: noninteracting bosons, the large U Bose-Hubbard model, 
and noninteracting fermions. We find that these parity correlations tend to be less strong than 
density-density correlations, but they carry similar information. 



I. INTRODUCTION 

Bakr et al. [l| recently reported site resolved images of 
atoms in optical lattices. Their tool, the quantum gas 
microscope, works by rapidly tuning the lattice near res- 
onance with an atomic transition: the resulting increase 
of the atomic confinement freezes the atoms in place. Si- 
multaneously, the sample is bathed with a separate near- 
resonant beam of light, causing the atoms to fluoresce, 
and allowing individual atoms to be observed. While 
this technique enables single-atom imaging it has a 
side effect that sites with an even number of atoms are 
quickly emptied by light-assisted collisions, and appear 
dark. The same effect causes all but one atom to be lost 
at sites with an odd number of atoms. Hence, the quan- 
tum gas microscope images map out the atomic number 
parity Px = (1— (— 1)"'')/2, where Ux is the density at 
site X. Here we investigate correlations 
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which can be measured by taking autocorrelations of the 
quantum gas microscope images. We compare this cor- 
relation function to the more typical density-density cor- 
relations, (nx7T,y) — (rix) (ny). Measurements of density 
correlations have been quite valuable for learning about 
cold gas systems and have future applications 

Of particular interest is a protocol introduced by Zhou 
and Ho[10] which allows one to accurately determine the 
temperature of a cold atomic system by averaging the 
density and density fluctuations of the system across the 
entire harmonic trapping region. 

We consider several different cases: (1) free or weakly 
interacting bosons, (2) the large U Bose-Hubbard model, 
and (3) free fermions. While D^y is much harder to 
interpret than the density-density correlator {nxTiy) — 
(rix) (riy), we find that it is a very useful probe. In par- 
ticular, there is much information contained in these cor- 
relation functions which is not simply found in the mean 
(Px). Although these are all simple models, our results 
are nontrivial. Calculating these correlations can be dif- 
ficult, even for noninteracting bosons. Furthermore, as 
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all of these systems have been realized, our results are 
relevant to experiments. 

There are many interesting systems (such as the t-J 
model) , where the occupation on a given site is restricted 
to be 1 or 0. For those systems there is a one-to-one 
correspondence between the parity correlations and the 
density-density correlations. This makes these models 
extremely well suited to experiments of this type, but we 
will not treat them here since their correlation functions 
can be calculated through standard means. 

Additionally, one may be able to use knowledge of the 
underlying physics to directly extract rii from the mea- 
surement Pi. For example, in the large U Hubbard model 
number fluctuations are sufficiently attenuated that at 
any point in space the density is likely to take on only 
one of two values. Depending on what quantities one 
is interested in, studying these images may be more 
fruitful than the analysis presented here. 

Throughout we will consider a uniform system. We 
envision using a local density approximation (LDA) to 
apply our results to a trapped gas. For the LDA to be 
valid we need two conditions. First the confinement of 
the cloud must be sufficiently weak that the system is 
locally homogeneous. Typically this condition is satisfied 
if the harmonic trapping frequency is small compared to 
the tunneling rate, Wc ^ t. Second, we require that 
the two points in our correlation function, x and y, are 
separated by a distance much smaller than the size of 
the cloud. Thus there is an upper limit on the length 
scale over which we can measure Dxy [12]. If -Dxy decays 
sufficiently rapidly, then we can completely characterize 
this function by experiments on a finite cloud. 

The remainder of this paper is organized as follows. In 
section [TT] we give results for weakly interacting bosons. 
In section IIIII we consider a bosonic Mott insulator. In 
section [TVl we discuss the two-component Fermi gas. 



II. BOSONS: IDEAL AND WEAKLY 
INTERACTING 

The most technically challenging of our calculations 
will be for the case of finite temperature noninteract- 
ing Bosons. We will also explain how weak interactions 
can be included, but will not give detailed expressions. 
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FIG. 1: Characteristic images of the density (left) and parity 
(right) of a homogeneous noninteracting BEC in an optical 
lattice at zero temperature with average density (n) = 2. 
Left: brighter colors correspond to higher density. Right: 
white sites correspond to odd numbers of particles, and black 
corresponds to even. 



In these systems the correlator D^y encodes information 
about the density-density correlation length and number 
squeezing. Our central result will be analytic expressions 
for D^y which are exact in the thermodynamic limit. 



A. Bose-Einstein Condensate in the ideal gas 

The simplest limit to consider is a zero-temperature 
ideal (noninteracting) Bose gas. As is well described in 
the literature ITS'], the number occupation on different 
sites will then be Poissonian, and uncorrelated, yielding 
D. 



0. The expectation value of the single site atom 



number parity operator is 
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)/2. 



(2) 



Figure [T] shows the expected quantum gas microscope 
images. 
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where Ck annihilates a boson with quasimomentum k. 
For calculating the properties of this system we will work 
on a finite lattice of N sites, and it will be understood 
that all calculations are to be evaluated in the limit 
N oo. Our results will apply to a wide varieties of 
lattices with a range of different parameters. We shall 
require that the dispersion ek is even, which is generi- 
cally true unless time reversal symmetry is broken. In 
this subsection we will further require that (ek — is 
positive definite, implying that the gas is noncondenscd. 
Under these conditions we will be able to calculate Z3xy 
as a function of /Lt and T. In section III B 31 we will analyze 
the superfluid state. 

We perform our calculation in a momentum occupation 
basis: 
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where M is the number of distinct occupied momentum 
states. The partition function can be written as 
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were the sum is over all possible occupations of the mo- 
mentum eigenstates. 



B. Finite temperature ideal Bose gas 

Scaling energies by the temperature, a gas of nonin- 
teracting bosons on an optical lattice is described by the 



1. Single Site Parity Expectation Value 



We begin by calculating the expectation value 
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where we define Ef^ = (e^ — fi) and /3 = 1/T. We define our fourier transform conventions as: 
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To evaluate the expectation value, we commute the operator (—1)"'' through all the operators to the vacuum state, 
where it evaluates to 1. In the position basis, (—1)"" = (1 — 2(5xy) (—1)"''. In the momentum basis, this yields 
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where A = 1 will be used as a formal expansion parameter. Using these operators we write 
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This expectation value can be calculated as a product of 
vacuum expectation values via Wick's theorem. We con- 
struct a formal power series in A [introduced in Eq. (|S])]: 
((—!)") = X]m^("^)'^™- construct this expansion 
by introducing the contractions 



Ck4 ^ [ck, dl] - <5kp = e'(k-p)-x. (10) 

Mathematically, F [m) counts the contribution to 
((—!)") in which there are m such contractions. Each 
of these contractions caries a factor of where N in 
the number of sites, but we will find an additional combi- 
natorial factor of 0{N"^) which offsets this. Thus, when 
A = 1 (as is physical), each term in this expansion will 
be of comparable size and we will need to sum the entire 
series to produce a useful result. 

We will begin by explicitly constructing F{m) for m = 
0, 1,2, and then will provide an argument for the general 
term. The leading term is F{Q) = 1, as the sum in Eq. ([9]) 
is then simply the partition function. 

Next, F(l) comes from a single contraction of equal 

momenta, Cpdp. Contractions which involve two distinct 
momenta such as Cpdj^. will leave unequal numbers of Cp 
and operators in Eq. ©. The expectation value of 
these will vanish unless an additional off-diagonal con- 
traction is made, contributing another power of A. Let 

p be the momentum involved in the contraction Cpdp. 
There will be rip such contractions, each of which con- 
tribute {—2X/N) to ((—1)""'). Summing over the mo- 
menta p, we see 

^(1) = -|E("p) = -2H- (11) 
p 

At second order, there are two classes of nonzero terms; 
we can either choose two contractions of equal momenta, 

= Cpdp X Ckrfjj. or we can contract two pairs of un- 
equal momenta in the form = Cpdj^ x c^dp. In both 
these cases, all the momentum dependent phases from 
the exponentials in ((TT|) cancel. This is a consequence of 
translational invariance, and will be true at each order. 
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FIG. 2: Diagrammatic representation of the expansion in A, 
with the third order term shown explicitly. The cluster de- 
composition principle applies once the terms are grouped into 
rings as described in the proof, and the combinatorics of the 
summations over momenta lead to the exponentiation of all 
nonvanishing diagrams. The numerical contribution of F (3) 
comes from summing over all momenta for each leg in the 
sum of ring diagram, taking into account the combinatorics 
described in the text. 

As long as p 7^ k, each of these two terms contribute the 
same amount: f^^ = = npnk{2X/N)'^ . The case 
where p = k can be ignored in the thermodynamic limit, 
as the sum of their contribution will scale as 1/7V. We 
sum over k and p, dividing by 2 to take care of over- 
counting, to arrive at: 

F(2) = 4(n)^ (12) 

Generalizing to the case of of m contractions is best 
performed by introducing a diagramatic language and 

reorganizing. A contraction Ckdp is represented by a ver- 
tex. Each vertex has two lines, representing the momenta 
k and p. Vertices which are connected share the same 
momentum. X"^F(rn) is given by a sum of all possible 
closed configurations of m vertices, attributing a factor 
—2X/N to each vertex, and a factor (rip) to each line. 
Any given diagram will have associated combinatoric fac- 
tors to eliminate double-counting. The result for F{1) is 
just a single loop. The two contributions to F{2), are the 
product of single loops, and the ring diagram with two 
vertices. The diagrams for F(3) are also shown in Fig. [21 
yielding F(3) = -8(n)^ 

In any term where there is a product of m identical dia- 
grams, there will be a combinatoric factor of (l/ml). This 
coefficient represents the fact that permuting the differ- 
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ent diagrams does not produce new terms. This prop- 
erty implies that the hnked cluster expansion applies, 
and upon exponentiating, disconnected diagrams factor 
out and can be discarded. Therefore, J2m ^"^F{m) — 
exp(E™^"^("^)), where R{m) = C™(-2A/A^)"(n)'", 
is the contribution from the linked diagram with m ver- 
tices: the TO- ring. The combinatoric factor Cm = 1/to 
results from the to different starting position of the ring. 
It can also be thought of as the ratio of the (to — 1)! ways 
of arranging m momenta into a ring, divided by the to! 
different possible permutations of those momenta. 

Since each TO-ring amplitude scales as the TOth power 
of the average density, our end result is a series in the 
density of the system. Summing the series yields 



((-1)-) = n 



1 



l + 2(n) 



In terms of our observable operators: 



Px = 



1 + 2 n 



(13) 



(14) 



When {n) <C 1, one has the expected behavior that 
('^x) = (-Px)- As (n) — ^ oo, there is no bias towards 
even or odd occupations and (Px) 1/2. 



2. Two-Site Expectation Value and Correlation Functions 

We now want to evaluate ((-1)"" (-1)"°)- AU of the 
combinatorical arguments of the previous section still 
hold, meaning that the expectation value can still be 
written as the exponential of the sum of ring ampli- 
tudes. However, beyond the 1-ring amplitude, the mo- 
mentum dependent phases will no longer cancel, making 
the ring amplitudes depend on x. Adopting the conven- 
tions above, we make the following replacements: 
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As in the previous section, A 1 is a for- 

mal expansion parameter. The diagram struc- 
ture is identical to the previous section, but each 
vertex now contributes (2A/iV)(l -I- e^^P"*")'^) = 
(4A/iV)e*(P-'')-="/2 cos ((k - p) • x/2). The ring diagrams 
can be written: 



Ri (x) 

i?2 (x) 
Rtu (x) 



-4(71), 

— ^ (nk^p) Cp_kCk- 



(16) 



k,p 



toAT" 



XI ("ki---nk„) 
ki-k„, 

xCk,„-ki ■ ■ ■ C'k„_i-k„ 



Cp = cos(p • x/2) 



We will integrate out one of the momenta in the expres- 
sion for Rm, generate a recursion relationship. In par- 
ticular, by using basic trigonometric identities we may 
write 



/kq - XI 



k+p 



N ----^-^-^ 2 2 
q 

where the Greens function Gx is defined by 

- E - (4-o) . (17) 



To produce a closed set of recursion relationships, we 
introduce a modified ring amplitude 



M^m(x) 



N 



2™ 

TO 



X ('^ki X ...nk,„) (18) 

ki...k„ 

Ck„_2-k„_i X Ck„_i+k„, 



where the final cosine depends on the sum of two neigh- 
boring momenta rather than the difference. If one repeats 
the same procedure for Eq. ((TS]), integrating out k^^i, 
one will find a term where two consecutive cosines depend 
on the sum of momenta. By taking k„i_2 ^ km-2, 
and using time reversal symmetry, one can express such 
a product in terms of Rm-i- The resulting recursion 
relations are 

Tn — 1 Tfi — 1 

Rr,r =-2(n) Rrn-1 - 2Gx W^-l, 



W„ 
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Tn - — 1 Tn — 1 

-2Gx i?™-i-2(n) Wm-i- (19) 



TO 



TO 



It is straightforward to solve these equations for Rm ± 
Wm, to find 

Rm±Wrn = (2/TO)[-2((n)±Gx)]'" (20) 

where we have used that i?i = — 4(ri) and Wi = — 4Gx. 
This yields the correlation function 



X^m = -log [il + 2{n)f -AGl 

m— 1 

((_!)»= (-1)"°) = ^ 



il + 2{n)r-AGl 
This allows us to express the correlator in Eq. ^ as 
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For the non-interacting gas, the density-density correla- 
tion function is simply the product of Greens functions. 
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larger than 1, implying that the parity correlations are 
always smaller than the density-density correlations. In 
particular, at large densities, Z^xy ^ {n) ^ ■ The charac- 
teristic length scale of these correlations, is however set 
by the density-density correlations. 



FIG. 3: Top: new vertex and line needed to include conden- 
sate. Bottom: Contribution to (( — 1)") = &q){^^Rm + Lm). 

Thus the parity correlation function is only a function 
of the average density (n) and the density-density corre- 
lation function. The denominator in Eq. p3p is always 



3. Base Condensed Phase 

The previous formalism can be modified slightly to de- 
scribe the Bose-Condensed phase, where there is macro- 
scopic occupation of the fc = mode. We work in an 
ensemble where the number of condensed particles uq is 
fixed and extensive, so 



((-!)-) = , 'm , (0| (^o)"" (ckj"^ • • • (-1)"- UT ■ ■ ■ Ur |0) e-^^ ^(25) 
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where all sums over k do not include k = 0. 

One readily includes the extensivity of uq through a 
modification of our diagramatic language. To motivate 
our formalism, we first consider the zero temperature ex- 
pectation value 
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In the thermodynamic limit (1— 2A/iV)"" — exp(— 2Ans), 
where the condensate density is Ug = uq/N. This agrees 
with our argument in Sec. Ill Al 

In the thermodynamic limit {nQ,N — >■ oo) this zero 
temperature gas has F{m) = (— 2A/A^)™no!/(m!(no — 
~^ (— 2Ans)™/m!. The natural way to express this 
diagrammatically is to again have m vertices, each of 
which contribute {—2X/N), and attach to each of them 
two short jagged lines, each of which contribute y^no. 
These jagged lines attach to only a single vertex, and in a 
coherent state formulation would have the meaning of the 
condensate order parameter (ao) f3|. The denominator 
ml represents the various permutations of the vertices. 

At non-zero temperature, one will have both contri- 
butions from contractions with non-zero k, described by 
the diagramatic language introduced in Sec. Ill Bl and 
contractions with k = 0. Generically one will have ring 
diagrams (Fig. [2]), and line diagrams, illustrated in Fig.|3l 
The diagrams describing Eq. (1^51) are a special case of the 
line diagrams. We again exponentiate the series, finding 
((-!)") = exp{J2m^m + Lm), where i?„ is the previ- 



ously introduced m-vertex ring-diagram (with the k = 
mode removed) and L,„ is the m- vertex line diagram. 



Lm — n-o 



-2A 



nki...nk„Cod]j^^ Ckidk2---Ck„rfS 



ki...k„ 

= Us (nns)" (-2A)" , 

where n^s — (n) — ng is the density of non-condensed 
particles. Summing this series yields: 



((-1)'^ 



1 



-2n„ 



1 + 2nn 



■ exp 



1 -I- 2nn 



(27) 



A similar approach can be used to calculate 
((—1)"" (—1)"°). As before, one simply changes the ver- 
tex to account for the extra phase factors. Our previous 
calculation of the ring diagrams goes through unchanged. 
For the line diagrams we again produce a recursion re- 
lationship by integrating out one of the momenta. The 
result is that 

Lm = ns (n„, -I- Gx)" (-2A)'" , (28) 

where Gx ~ Gx — ng. The correlator then becomes 



((-ir (-ID 



exp 



-4ns 



(l+2(n„, + G.c)) 

(1 + 2n„,)2 - 4G2 



(29) 



The even-odd correlation function will thus decay expo- 
nentially to zero with increasing superfluid density. 



4- Weak Interactions 

The diagramatic series which we have developed may 
readily be expanded to allow for a perturbative treatment 
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of interactions. A new vertex, representing interparticle 
interactions will appear. Physically, repulsive interac- 
tions will suppress number fluctuations, hence increase 
((— 1)""). One expects that when temperature is large 
compared to the interaction energy the correlation length 
will still be dominated by thermal effects, and our free 
particle results will be valid. The lowest temperatures 
produced in experiments are a fraction of the bandwidth, 
and hence a fraction of the hopping t. Thus when t^U 
(where U is the on-site interaction energy) our theory 
should suffice for describing the correlations. Moreover, 
a Feshbach resonance can always be used to reduce the 
importance of interactions. 



III. LARGE U BOSE-HUBBARD MODEL 

A. Perturbation Theory about the Mott state 

We now consider the opposite limit, where interactions 
are strong compared to the hopping. This is typically 
described by the Bose-Hubbard model [13,] : 



Hbh = 2^ {Urii {n, - 1) - fMUi) 



Here, fii is a site dependent potential that combines 
the bare chemical potential and any trapping potentials. 
Here we will just treat the uniform system, taking /i^ = /i 
- as before we will use a local density approximation to 
treat the trap. We want to calculate correlation functions 
of the form: 

Dju ^ {P,Pk) - (Pj) (Pk) , 

where Pj is the even-odd projector. 

Following Freericks et aZ[l5|, and the related works 
of Eckardt et aZ[l^, we perturbatively calculate the 
nearest-neighbor correlator in powers of t/U. This ex- 
pansion breaks down in the superfluid phase, but cap- 
tures the leading order correlations in the Mott phases. 
The other typical approximation used to investigate the 
Bose-Hubbard model, the Gutzwiller ansatz, is not ap- 
propriate for calculating Djk- By fiat, the Gutzwiller 
ansatz restricts correlations to be zero or infinite range. 
We will only work to second order in t, although the 
generalization to higher orders is straightforward. 



Second order perturbation theory gives 



(P,) - (P,-)o = f dT [ dT'{H^op{T)H^op{r')P,)Ho 

JQ Jo 

dT'{Hhop{T)Hhop{T'))Ho, (31) 
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dT 

^0 



where Hq = J2i (U n-i (n. 

H, and Hhop = - Y^a ^ 
the local ensemble is 



- 1) — ^rii) is the local part of 
^ " ■ The expectation value in 



(32) 



Tre^/^ffo ' 

and the interaction picture operators Hhop[T) = 
e~^oT jj^^^^Hqt ^ p^-^ expression similar to (1511) describes 

{PiPk)- 

In Eq. (PT|) . we introduce a resolution of the identity 
between each operator, consisting of the states with a 
fixed number of particles on every site. As detailed in 
[isj the resulting expression factors, leaving only chains 
of hopping operators which pass through site j (or k in 
the case of calculation (PjPk)). The resulting expressions 
are most conveniently written in terms of the local Greens 
functions 

G,{T,T')^-{T^a,{T)a]{T')) , (33) 

\ / Ho 

Gf (r, r') ^ - (Tra, (r) a] (r') pA (34) 

\ / Ho 

where the time ordering operator Tr places operators at 
earlier imaginary times towards the right. Implicitly we 
set the imaginary time of the P operators to zero so that 
they always appear on the right. We define the bare 
single site energies, e„ = (Un (n — 1) — /in), and the sin- 
gle site partition function Zgs — X^n^"'^'"- The Greens 
functions are then explicitly given as 



G3 (r, r') ^Y.P- [(" + 1) e (r - r') ei^' 

n 

+ Y.P^ p(r'-T)e(^-^')^.- 



(35) 



where p„ = e ^"^/Zss, and = e„±i - e„. {t,t') 
is simply given by removing the even n terms in the ex- 
pression for Gj . One then finds 
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(P,) - {Pi), 



vPP 



np 1 ^ 

dx dy Gk {x, y) Gj [y, x) = 2_, P"/'™ (^nm + ^mn) , 

-'0 ^0 ^ ^ 



/ da; / dyG^ {x,y)Gj{y,x) ^y^PnPmP {n){Tnrn+'^nin) , 

I dx I dy Gf (x, y) Gf (y, x) = ^ PnPmP (") (m) (r„„i + r„i„) 
(1 + to) n (-1 + /?[/ (1 - m + n)+ e^C^(-i-™+")) 

(1 + TO — ?l)^ 



As shown, the integrals reduce to a rapidly converging 
double sum which we compute numerically. We plot 
these correlation functions in Fig. |4l The dominant fea- 
ture is a series of plateaus where Djk is positive, punctu- 
ated by downward spikes where Djk becomes negative. 
The size of the plateaus increase monotonically as /i in- 
creases. 

Both the plateaus and the spikes can be qualitatively 
understood from a truncated two-site model. Deep in 
the TO-particle Mott regime, one expects that config- 
urations at two neighboring sites would be dominated 
by the states {|to, to) , |to — 1, to -I- 1) , |to -I- 1, to — 1)}, 
where the two integers represent the number of parti- 
cles on each site. The probability of being in either of 
the latter two states is 6 ^ m(t/U)'^, where the factor 
of TO comes from Bose enhancement. The correlations 
function (PjPk) - {Pj){Pk) = 4(6- b^) is positive, re- 
flecting the fact that the parity of two neighboring sites 
is correlated. 

Conversely, near the single-site level crossings {fJ./U — 
to), the configurations of neighboring sites are dominated 
by {|to, to) , |to, to -|- 1) , |to -|- 1, to)}. At finite t the lat- 
ter two states are energetically favorable, resulting in an- 
ticorrelations between the parity at the two sites. Al- 
though in this regime our perturbation theory breaks 
down, the anticorrelations should be robust. 



Djt(t=0.03U, T=t,2t,4t,8t) 




FIG. 4: (Color Online) Leading order strong coupling per- 
turbation theory calculation of nearest neighbor parity cor- 
relations Djk at temperatures T = {t, 2t, it, 8t} (from top to 
bottom in the plateaux; blue, purple, yellow, green in color 
versions of this article) for t = 0.03C/ as a function of fJ^/U. At 
small t/U these are positive in the Mott state, and negative in 
the superfluid. The perturbation theory breaks down in the 
superfluid, but the sign of the correlations should be robust. 



and the parity correlations can be readily expressed in 
terms of density correlations. In the noninteracting limit, 
these factor, yielding for x 7^ 



IV. SPIN-i FERMIONS 

We readily calculate the parity correlations of spin- 1/2 
fermions. Since each site can have at most one particle 
of each spin, the parity is 

Pi = + nil - 2"itn4, (36) 



(P.Po) - (Px) (Po) = 4G4 - 2 (1 - n/2f G^ , (37) 

where G^ = ^ J2k ("k) e"^^'^ is the single particle den- 
sity matrix. The correlations vanish when n = 0, 2. For 
intermediate n the correlations can be positive or nega- 
tive at short distances, but are always negative at large 
distances. 
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V. CONCLUSION 

We have calculated the atomic parity correlation func- 
tions D^y for three of the most common systems in op- 
tical lattice physics, non-interacting bosons, the Bose- 
Hubbard model, and non-interacting fermions. In the 
free boson system, we found that Z^xy decays to zero 
at long distances at the same rate as the square of the 
thermal Green's function, matching the behaviour of the 
density-density correlation function. In addition, the 
magnitude of D^y decays to zero with increasing density, 
and with increasing condensate density. In the Bose- 
Hubbard model, we demonstrated a simple link between 
the parity and density operators, and calculated the par- 
ity correlation functions to second order in t/U using 
imaginary time perturbation theory. 



In fermionic systems, the even-odd correlators are cal- 
culated in the same manner as the density-density cor- 
relation functions, and no new theoretical machinery is 
needed to interpret them. 
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